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Abstract. Several instrumental effects of the Long Wave- 
length channel of ISOCAM, the camera on board the In- 
frared Space Observatory, degrade the processed images. 
We present new data processing techniques that correct 
these defects, taking advantage of the fact that a position 
on the sky has been observed by several pixels at different 
times. We use this redundant information (1) to correct 
the long term variation of the detector response, (2) to 
correct memory effects after glitches and point sources, 
and (3) to refine the deglitching process. As an example 
we have applied our processing on the gamma-ray burst 
observation GRB 970402. Our new data processing tech- 
niques allow the detection of faint extended emission with 
contrast smaller than 1% of the zodiacal background. The 
data reduction corrects instrumental effects to the point 
where the noise in the final map is dominated by the read- 
out and the photon noises. All raster ISOCAM observa- 
tions can benefit from the data processing described here. 
This includes mapping of solar system extended objects 
(comet dust trails), nearby clouds and star forming re- 
gions, images from diffuse emission in the Galactic plane 
and external galaxies. These techniques could also be ap- 
plied to other raster type observations (e.g. ISOPHOT). 

Key words: Techniques: image processing - Infrared: gen- 
eral - Instrumentation: detectors 



1. Introduction 

The Infrared Space Observatory (ISO) of the European 
Space Agency died on April 8 1998, 10 months after its ex- 
pected end. This mission has been a great success (Kessler 
1999) and many of the scientific discoveries are yet to 
come. Much of the scientifical analysis of the data are 



Send offprint requests to: Marc-Antoine Miville-Deschenes 

* Based on observations with ISO, an ESA project with in- 
struments funded by ESA Member States (especially the PI 
countries: France, Germanv. the Netherlands and the United 



still held by data reduction problems. In particular for 
the camera on board of ISO (ISOCAM - Cesarsky et al. 
(1996)), instrumental effects which are not well under- 
stood limit the sensitivity. 

Much effort has been made to model the response of 
the ISOCAM array based on theoretical understanding 
of infrared detectors (Abergel et al. (1999), Coulais & 
Abergel (2000) and references therein). Nevertheless, no 
overall model of the ISOCAM response to flux steps and 
glitches is available. ISOCAM data reduction thus relies 
on an empirical understanding of the detector. Recently, 
sophisticated data reduction techniques have been devel- 
opped to take into account some aspects of ISOCAM re- 
sponse (Starck et al. (1999), Desert et al. (1999), Aussel 
et al. (1999), Altieri et al. (1999)). These methods are 
close to being optimal for the detection of point sources. 
Nevertheless, in many observations, instrumental effects 
are still preventing the study of faint extended emission. 

The development of the data reduction method pre- 
sented in this paper was motivated by the analysis of 
ISOCAM observations of diffuse and translucent inter- 
stellar clouds. Such clouds, when illuminated by the so- 
lar neighborhood radiation field have mid-infrared bright- 
ness with very low contrast with respect to the zodiacal 
light (at most a few percent). To unveil these low signal- 
to-noise ratio clouds, variations in the detector response 
have to be corrected to an accuracy better than a fraction 
of 1%. At the present time, data processing techniques 
available do not allow to achieve such high sensitivity for 
extended emission. In this paper we present a data pro- 
cessing method that make it possible. 

The algorithms presented here apply to observations 
made in raster mode with the Long Wavelength channel 
of ISOCAM (LW). It makes use of the fact that a given po- 
sition on the sky has been observed several times by differ- 
ent camera pixels. Tests have been performed on extended 
emission observations presenting low and high contrasts. 
We have checked that our methods can be applied to most 
raster mode observations, and as a consequence concerns a 
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To illustrate the data processing we use two different 
observations of the same field, obtained subsequently in 
exactly the same configuration. The amplitude of instru- 
mental effects are not the same in both observations, giv- 
ing us constraints on the validity of our method. In section 
§ [2| we present the observations used to illustrate the data 
reduction chain. In section § [| standard data reduction 
techniques are briefly presented. The main problems of 
the standard reduction are described in § ^ and the new 
data processing approach to address these problems is de- 
tailed in § ||. The performances of the overall method are 
discussed in § ||. 

2. Observation 

The main ISOCAM technical characteristics are presented 
in Cesarsky et al. (1996). The LW channel of ISOCAM 
operates between 4 and 18 ^m. A lens wheel allows the 
selection of the field of view per pixel (1.5, 3, 6 and 12 arc- 
sees) , and a filter wheel the selection of the spectral band 
pass (10 broad band filters and two Continuously Variable 
Filters). The detector is a 500 /im thick crystal made of 
Gallium doped Silicon photo-conductor hybridized by In- 
dium bumps. The 32x32 square pixels are defined with a 
pitch of 100 /im. 

A given observation (characterized by a given config- 
uration (lens, filter, integration time)) is presented as a 
collection of 32 x 32 images (called readouts) gathered to- 
gether in a data cube. Three different observational modes 
were available with ISOCAM (Siebenmorgen et al. 1997): 
1) single pointing, raster mode, 2) beam switching and 
3) spectrophotometry. The observations presented in this 
paper were obtained in raster mode where many readouts 
are put together to build a mosaic of the observed object 
(called sky image). The scanning strategy is made in such 
a way that each readout taken on the sky has some overlap 
with its neighbours (see Figure 0) . 

To illustrate the data reduction procedures, we use 
two ISOCAM observations of the gamma-ray burst GRB 
970402 first observed by BeppoSAX (Feroci et al. 1997) on 
April 2.93 UT 1997. The coordinates (J2000) of the burst 
were a = 14 /l 50 m 16 s , 6 = -69°19'.9 with an error circle 
of radius 3'. Target-of-opportunity ISO observations (ISO- 
CAM and ISOPHOT) were requested by Castro-Tirado 
et al. (1998) to detect a transient infrared emission of 
this burst. Unfortunately the observations did not show 
such emission. Nevertheless these particular observations 
are of great interest for us since the same field has been 
observed twice, subsequently within the same orbit, in ex- 
actly the same configuration to look for a decrease in the 
GRB emission. The observational strategy of the obser- 
vations is sketched in Figure [l] and described in Table |l|. 
Comparing the result of our processing on both obser- 
vations (GRB1 and GRB 2 in the following) allows us to 
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Fig. 1. Raster observing strategy. The crosses indicate 
the positions on the sky of the camera center. The first 
(solid line) and the sixteenth (dash line) positions of the 
camera are shown as square. The arrow shows the scanning 
direction. 



3. Building the standard low level sky map 

The standard ISOCAM reduction method consists of the 
following steps: 

Deglitching. On any observation, numerous energetic 
particles impacts leave traces on one or several pixels of 
the detector. Most of these impacts do not affect the de- 
tector response on long time scale and are seen as in- 
stantaneous flux step (fast glitches) . In Figure |2[ a typ- 
ical flux history of one pixel as function of time is shown 
where many fast glitches are apparent. Various techniques 
achieve to suppress these emission spikes (Starck et al. 
1999). 

Dark image subtraction. The response of the LW 
detector of ISOCAM is not zero when the detector does 
not receive any photon coming from the sky. This dark im- 
age must be subtracted from the data. Temporal variation 
of the dark image has been observed (Starck et al. (1999)) 
and modeled by Biviano et al. (1998) using a polynomial 
function of time. 

Short term transient. It has been known since pre- 
launch measurements (Perault et al. 1994) that the ISO- 
CAM detector is affected by a transient effect, i.e. its re- 
sponse to a given flux step is not instantaneous. Recently 
a new model of the response, based on the physic of Si:Ga 
arrays (Fouks & Schubert 1995), has been developped by 
Coulais & Abergel (2000). The uncertainty on the cor- 
rected flux with this new model is ~ 1% but is raises 
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GRB1 


GRB2 


Date 


April 5 1997 


April 5 1997 


ISO Revolution 


506 


506 


Time since activation 


2.12 (hrs) 


3.42 (hrs) 


Filter 


LW10 


LW10 


Wavelength 


8.5-15.0/im 


8.5-15.0/im 


Right Ascension J2000 (center) 


14h 50m 18.2s 


14h 50m 18.1s 


Declination J2000 (center) 


-69° 19' 57.1" 


-69° 19' 59.3" 


Total number of readouts 


907 


907 


Number of raster steps 


8x8 


8x8 


Readouts per positions 


12 


12 


Step size (pixels) 


8 


8 


Exposure time (seconds) 


5.04 


5.04 


Pixel size 


6" 


6" 


Comment 


Long term transient 


Good quality 



1000 2000 3000 4000 

t (sec) 

Fig. 2. Temporal evolution of a typical pixel. A) Raw 
data - many energetic electrons have hit the pixel causing 
instantaneous flux steps known as fast glitches. A response 
accident is apparent near t = 350 seconds caused by an 
ion impact (slow glitch). B) Fast glitches corrected data - 
short time flux steps have been identified as glitches and 
removed from the data cube. 



strong problem here since we are mostly interested in dif- 
fuse emission. 

Flat Field. Generally, each pixel of a detector does 
not have exactly the same response to a given flux. This 
spatial variation of the detector response (the flat field 
in the following) must be taken into account. The most 
straightforward way to compute an image of the flat field 
is to average all readouts of a data cube, and normalize the 
mean to 1 (Starck et al. (1999)). This technique supposes 
that all camera pixels have seen the same average flux 



the sky much larger than the field of view of the camera, 
which do not contain any systematic gradient. 

Building the sky image. After dividing the data 
cube by the flat field, the final standard data reduction 
step is to project each readout on the sky to build the 
sky image. The LW channel of ISOCAM is affected by a 
field of view distortion problem that must be taken into 
account in this operation (Starck et al. 1999), especially 
to co-add properly emission from point sources and small 
scale structures. 

The sky image of the GRB1 observation obtained with 
the standard data processing described in this section is 
shown in Figure |^a. It is affected by many instrumental 
problems that prevent to take advantage of the full sensi- 
tivity of the instrument. In the next section we describe 
the open problems for the imaging of faint extended emis- 
sion. 

4. Open problems 

Long term transient. Figure || presents the temporal 
evolution, of the median flux observed in the 11x11 cen- 
tral square of the detector for the GRB1 (a) and GRB2 (b) 
observations. One striking feature of Figure || is the drift 
in the first observation, not seen in the second one. This 
drift in ISOCAM response is observed quite frequently and 
is called in the following the Long Term Transient (LTT) . 
When observed, the amplitude of the LTT is about 5% of 
the total sky emission like in Figure One also sees that 
the first sky image (Figure ||a) is completely dominated 
by the long term drift. This is true for any observation af- 
fected by a LTT of faint extended emission with contrasts 
smaller than a few percent of the zodiacal light. 

The LTT has been identified during the pre-launch 
measurements. Many calibration observations present a 
slow and continuous increase of the response after a posi- 
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Fig. 3. Temporal evolution after short term transient cor- 
rection of the median flux observed in the 11x11 central 
part of the camera, for the GRB1 (A) and GRB2 (B) ob- 
servations. 

the short term transient, no analytical or physical descrip- 
tion has been developped up to now to correct this spe- 
cific effect which seems to be a general behavior of this 
type of detector. Vinokurov & Fouks (1991) adressed the 
problem of the long term drift in their study of nonlinear 
response of photoconductor. Unfortunatey the analytical 
approximation they propose does not seem to reproduce 
accurately ISOCAM's LTT. 

Slow glitches. Another important problem is related 
to glitches which modify the response during a significant 
time. This is illustrated in Figure || where one notices a 
significant change in the detector response near t = 350 
seconds, just after a glitch impact. The glitch impact is 
removed by the standard deglitching algorithm but the 
detector response is significantly disrupted for more than 
500 seconds. In this example the signal is depressed after 
the glitch, but it can also be enhanced. These glitches with 
memory effect (called in the following 11 slow glitches"), 
which affect the response of one or more pixels for some 
time after the impact, are believed to be due to heavy ions. 
They are responsible for most of the periodic patterns seen 
in Figure ||a. No model has been developped to correct this 
instrumental effect. Desert et al. (1999) describe the vari- 
ous types of slow glitches known and a method to correct 
the response accident. However, this method has been op- 
timised for point sources extraction project where diffuse 
emission is not intented to be restored, so we do not use 
it. 

Ghosts. There are also several artificial point sources 
("ghosts" in the following) due to uncorrected memory 
effects. After seing a bright point source, the response of 
a given pixel is significantly affected for some time. We 
have seen in the previous section that significant memory 



moves on the sky from one sky position to the other, pixels 
that have observed bright point sources are affected by 
a memory effect, and the pattern of each point source 
appears repetitively in the sky image until the response of 
the pixels gets back to a "normal" value. 

We conclude that the actual data processing does not 
adequately take into account (1) the long term drift, 
(2) slow glitches and (3) short term transients on point 
sources. The long term drift precludes the study of large 
scale emission fainter than 10% of the zodiacal light back- 
ground. All instrumental effects limit the brightness sen- 
sitivity for small scale structures. Solving these problems 
is crucial to take advantage of the unique sensitivity and 
angular resolution of ISOCAM. This is particularly true 
for small scale structures since generally, the contrast of 
interstellar medium emission decreases from large to small 
scales (Gautier et al. 1992). 

5. New approach for LW-ISOCAM data 
processing 

5.1. General description 

5.1.1. General inversion method 

To address the problems described in § [| we have devel- 
opped a method which uses the fact that a position on 
the sky has been observed by several ISOCAM pixels at 
different times (see Figure |l|) . The redundant information 
is used to separate the contributions of the sky emission 
and of instrumental effects to the observed signal. This ap- 
proach has already been used for the IRAS mission (e.g. 
Okumura (1991) and Wheelock et al. (1997)) and may be 
generalized to every raster type observations, whatever the 
wavelength of observation. Formally the processing of as- 
tronomical data with spatial redundancy could be treated 
as an inversion problem. We can consider that the data ob- 
served O is the result of the convolution of the real sky S 
with the instrumental function / plus some additive noise 
N: 



= I<E)S + N. 



(1) 



When the observation has been conducted with spatial 
redundancy, it may be possible to address the data pro- 
cessing problem globally by finding /, S and N that will 
reproduce O and minimize the difference between mea- 
surements obtained at the same sky position. Such an ap- 
proach supposes that the main instrumental effects that 
affect / and N are known enough to constrain the inver- 
sion method. Presently this is not the case for ISOCAM. 

5.1.2. Sequential approach 

ISOCAM's response variations are complex and presently 
we are not able to model it with a reasonable number 
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READ DATA 



DEGLITCHING 



(these are called bad pixels in the following) . This last op- 
eration removes slow glitches and ghosts. These two op- 
erations (variable flat field and bad pixels identification) 
can be done several times to improve the sky image (see 
Figure ||) . 



DARK SUBSTRACTION 



SHORT TRANSIENT CORRECTION 



LONG TRANSIENT CORRECTION 



BUILD SKY IMAGE 



TIME DEPENDANT FLAT FIELD 



BUILD SKY IMAGE 



BAD PIXELS IDENTIFICATION 



Fig. 4. Data reduction chain illustrating the data pro- 
cessing steps from top to bottom. The time dependant 
flat field and bad pixels identification operations can be 
done iteratively. 



adopted a sequential approach where the instrumental ef- 
fects are treated one at a time. Nevertheless, even if we 
cannot use such a global method, the idea of minimizing 
the differences between pixels that have seen the same sky 
position, is the milestone of our approach. In particular the 
ISOCAM long term drift problem is addressed by a least 
square minimization technique based on the fact that, if 
the detector is slowly reaching stabilization, the measured 
flux is also approaching the real observed flux. Here we 
suppose that every pixels of the array is affected by the 
same long term drift. The hardness of this assumption will 
be discussed in § ||. 

Once the LTT is removed, the other response varia- 
tions (slow glitches, ghosts and residual transient effects) 
are corrected by comparing the data cube to the sky im- 
age. This is done by two operations. First we compute 
a variable flat field that takes into account pixel-to-pixel 



5.2. Long term transient 
5.2.1. Gain or Offset effect ? 

First of all it is important to determine how the long term 
transient, the flat field F and the observed flux I obs are 
related to the incident flux I s ky We have considered the 
three following possibilities for the LTT: 



1. Gain effect I obs = G(t) F I sky 

2. Offset effect affected by flat: I obs = F (I sky 

3. Pure offset effect: L obs — F I s k y + A(t) 



A(t)) 



The two GRB observations are of great help to con- 
clude on the form of the LTT. As the two observations 
were done exactly in the same way, it is possible to sub- 
tract directly the two data cubes. In Figure ||a we show 
the flux history of two pixels of the GRB2 observation. 
As one can see there is no long term drift detectable in 
this observation. The flux difference between the two se- 
lected pixels is due to a ^25% flat field difference. The 
difference between observation I and 2 for these two pix- 
els is shown in Figure ^d. Both pixels have a very similar 
drift; the 25% difference between the two pixels flux does 
not appear in figure |B|b. We then conclude that the LTT 
does not depend on the signal FI s k y . In this context, the 
only valid description of the LTT is the pure offset effect: 
Iobs — F I s k y + A(t). In the following we consider the LTT 
as a single offset over all the detector; in other words we 
suppose that all pixels are affected by the same offset. 

5.2.2. Formalism 

For a given pixel at a position (x, y) on the detector array 
and at a given time t, the observed flux I t, s (x,y,t) is 
related to the temporally varying flat field F(x,y,t), the 
incident flux I s k y (x,y,t) and the long term drift A(t) by 
the following equation: 



I obs (x,y,t) = F(x,y,t)I sky (x,y,t) + A(t). 



(2) 



Here we suppose that A(t) is not pixel dependant. 

The offset function A(t) is found using equation^ and 
the spatial redundancy inherent to raster mode observa- 
tions. We determine A(t) by solving a set of linear equa- 
tions obtained by comparing flat field corrected intensities 
of the same sky positions but obtained at different times. 
The flat field F(x, y, t) is computed using the perturbated 
single flat field method (see 
intensities can be written: 



A.2). The flat field corrected 
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Fig. 5. (A) Flux history of two pixels of the GRB2 obser- 
vation. This observation is not affected by the LTT. The 
flux difference between the two pixel is due to an intrin- 
sic response difference (flat field). (B) Difference between 
the two GRB observations for the two pixels considered in 
(A). Even if the two pixels have very different response, 
the LTT amplitude is the same. 



In raster mode the camera does not always point at 
the same position on the sky. To compare two pixels that 
have seen the same position on the sky we must put all 
readouts in a common spatial reference frame. Practically 
we project I obs (x,y,t) and F(x,y,t) on the plane of the 
sky (right ascension (a) - declination (5) reference frame) . 
Following equation ^, the difference between two pixels in 
that reference frame that have seen the same sky I s k y (a, S) 
at different time (ti and tj) is: 



I obs (a,S,U) I obs (a,5,tj) _ A(U) 



A(i, 



F(a,S,U) F(a,6,tj) F(a,8,U) F(a,5,tj) 



(4) 



To ensure that we compare fluxes obtained at the same sky 
position, the field of view distortion must be accurately 
taken into account in the projection operation. 

The function of interest A(t) is estimated using a least- 
square minimization technique with the following mini- 
mization criterion: 



E 



a,6,U,t 



I obs (a,S, U) _ I obs (a,S, tj) 
F(a,S,ti) F(a,S,tj) 

A(U) A(/,j) 



F(a,S,U) F(a,S,tj) 



(5) 



Here the sum is on all the possible pixel pairs that have 
seen the same sky. 

The function A(t) which minimizes the value of \ 2 is 
found by solving the system determined by: 



Equation ^ represents a standard set of linear equa- 
tions which can be written in a matrix form: AA(t) = B, 
with: 



Mh3)i& - E F(a,S,U)F(a,S,tj) 



A(i,j) i=j = 



*«= E 



ff F(a, 5,U)F(a,5,tj) 

Lbs {a, 5, tj) I obs (a,S,ti) 



F(a,6, U 



F(a,6,tj) F{a,5,ti) 



(7) 



(8) 



.(9) 



The A(i,j) and B(i) terms are computed with all pix- 
els that overlap. For each pair (i,j), the sums are always 
computed on the positions (a, 6) where I obs (a, 6, ti) and 
Iobs(ct, 5,tj) are both defined. 

Finally A(t) is found from: A(t) = A~ X B. As the sec- 
ond derivative of x 2 is always positive, the solution found 
for A(t) necessarily minimizes the % 2 criterion. We add an 
offset to all values of A(t) to force the correction to be zero 
at the end of the observation: A(t en d) = 0- This is justified 
by the fact that the long term transient tends to stabilize 
at the end of an observation. However, there are several 
observations not stabilised at the end, so that the absolute 
brightness may be systematically shifted (generally below 
a few %). For the observations of low contrasted clouds 
on top of a flat large scale emission (at least the zodiacal 
emission), this point is not critical since we always remove 
the large scale emission at the end of our processing. 

5.2.3. Practical implementation 

Beside detector noise, two additional sources of uncertain- 
ties affect the comparison of raster images and thus the 
LTT correction through the minimization algorithm de- 
scribed in the previous section: slow glitches and flat field 
variations along the observation. The signal measured on 
point sources is not identical in the individual readouts, es- 
sentially due to the undersampling of the the point spread 
function for the 6" pixels. Bright sources are thus an ad- 
ditional source of error but they can easily be discarded. 
The two other noise sources make the practical implemen- 
tation of the formalism a non-straightforward procedure. 

An extensive use of the LTT inversion method pre- 
sented here on real data has demonstrated the extreme 
importance of an accurate flat field. As the LTT correction 
is based on the comparison of the brightness measured by 
different parts of the array, its result depends on the accu- 
racy of the flat field. In particular, when F(a, 6, t) is not 
well estimated, oscillations in phase with the rastering of 
the observations may appear in the correction found. To 
pass around this difficulty we compute an approximate 
LTT correction (as oppose to the exact correction which 
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in some cases, the LTT can be approximated by the sum 
of two exponential functions, one upward and one down- 
ward: 



A(t) w P x exp (-Q xf fl )-Sx exp (-T x < c 



(10) 



where P, Q, R, S, T and U are strictly positive. By ap- 
proximating the LTT in that way, we greatly simplify the 
problem as A(i) depends now only on six parameters. As 
for the exact solution, the approximated solution is found 
by minimizing the criterion of equation ^|, using an IDL 
curve fitting program (MPFIT). The obvious advantage of 
this method is that the approximated LTT correction is 
smooth, getting rid of the oscillations found in some exact 
solutions. 

On the other hand, we have observed cases 
where this approximation of the LTT does not 
hold. Sometimes, the LTT shows some oscillations 
that can dominate the emission in low contrasted 
fields (e.g. cirrus clouds - see Miville-Deschenes 
et al. (2000)). In these cases, an accurate estima- 
tion of the flat field and the use of the exact solu- 
tion for the LTT are mandatory. 

5.2.4. LTT correction for the GRB1 observation 

Figure ^| presents the LTT corrections found for t he G RB1 



A. 2 ) was 



observation. A perturbated single flat field (see 
used and point sources were discarded to compute the ex- 
act and approximated corrections. The approximated LTT 
correction is smoother than the exact one, which oscillates 
with a ~ 0.2 ADU/g/s amplitude and a period of ~ 1800 
seconds. In this case we have applied the approxi- 
mated LTT correction. The sky image computed after 
that correction has been done is presented in Figure |^b. 
One sees that it is mandatory to apply the LTT correction 
since its amplitude (3 ADU/g/s) is nearly three times the 
amplitude of the emission of interest (1.1 ADU/g/s). At 
this stage we have used a single flat field to compute the 
sky image. It is necessary to use a variable flat field to 
correct the artefacts (e.g. periodic patterns) seen in Fig- 
ure ||b. 

5.3. Variable Flat Field 

Now that the LTT has been corrected, we take into ac- 
count the pixel-to-pixel temporal variations of the detec- 
tor response. These response variations are observed at 
various timescales. At short timescales they are due to 
bad short term transient correction, particularly on point 
sources. On longer time scales they are mainly caused by 
slow glitches (see § ||). The use of a single flat field (see § ||) 
does not take into account these temporal variations (that 
represent ~l-3 % of the average flat field) preventing to 
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Fig. 6. Exact and approximated LTT correction for the 
GRB1 observation. Estimating the error on the LTT cor- 
rection is difficult as it is very sensitive to the accuracy of 
the flat field used. The error bar represented here is the 
statistical error (0.08 ADU/g/s - see g). A more realistic 
error determination would be related to the accuracy of 
the flat field. For the GRB observations, a 1% uncertainty 
on the flat field corresponds to a 0.4 ADU/g/s error on 
the LTT correction. 

pixel response variations with a time-dependant flat-field 
(called "variable flat field 1 in the following). 

For LTT corrected data, the observed flux I b s at po- 
sition (x, y) on the array and at time t is 



I obs (x,y,t) = F(x,y,t) I sky (x,y,t), 



(11) 



where F(x,y,t) is the instantaneous response of pixel 
(x, y) at time t. Thus flat field and sky structures are 
mixed together in I f, s (x,y,t). In the following we show 
how the flat field variations can be extracted from the 
data by estimating I s k v {x,y,t) and by taking advantage 
of the spatial redundancy. 

In raster observation mode, many pixels of the data 
cube have seen the same position (a, S) on the sky. By 
averaging all these pixels we reduce the noise due to 
instrumental effects on the computation of the sky im- 
age I a ky (a, 6). The estimate of I s ky(x,y,t) is made by 
an inverse projection of I s k y (ot,5) on each readout of 
the data cube. Then P(x, y, t) is computed by averaging 
lobs (x, y, t) 1 1 s k v (x,y,t) with a sliding window on the time 
axis (see § |A.l| ). Here are the guidelines of this method: 

1. Construct a sky image. 

2. Smooth (median smoothing) the sky image with a 10 x 
10 window. 

3. Compute an ideal cube I s k y (x,y,t) by projecting the 
smoothted sky image on each readout of the data cube. 

4. Smooth (median smoothing) I t> s {x,y,t)/I s k y (x,y,t) 
on the time axis. The size of the smoothing window 
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Fig. 7. Ghost identification. Flux history of a pixel which 
has observed a point source. The flux values after the point 
source are affected by a memory effect and are identified 
as a ghost (flagged out). 



The sky image of the GRB1 observation, obtained with 
the variable flat field, is shown in Figure ||c. The size of 
the filtering window on the time axis is 100 (correspond- 
ing to ^500 seconds). As seen in Figure ||c, the variable 
flat field removes almost all periodic patterns due to high 
frequency variations of the detector response. Compare 
to other methods we have tested (see A.l and A.2), this 
variable flat field gives by far the best results. 



5.4- Bad pixels identification 

The deglitching process used at the beginning of the re- 
duction chain (see § ||) removes extremely deviant pixels 
that have been hit by cosmic rays. To go further in the 
noise minimization process, we take again advantage of 
the spatial redundancy. 



5.4.1. Ghost removal 

Memory effects often appear after a strong flux step (e.g 
after a point source) due to improper short term transient 
correction; this is what we call ghosts. To identify pixels 
affected by such effects, we look on the flux history of every 
pixel for memory effect after a flux step (see Figure Q). 
Again we use the redundancy information to improve the 
identification of ghosts. Here is how we proceed: 

1. Smooth (median smoothing) the sky image with a 3 x 
3 window. 

2. Compute an ideal cube by projecting the smoothted 
sky image on each readout of the data cube. 

3. Compute the residual cube = Idata cube - ideal cubel. 



5. On the time history of each pixels, identify ghost as 
residual fluxes above or under the noise level after a 
strong flux steps. 

6. Flag ghosts in the original data cube. 

5.4.2. Noise reducer 

We can go further in the reduction of the noise level by 
working sky position by sky position instead of working on 
the time history of every pixel. The idea is to look at pixels 
in the data cube that have seen the same sky position and 
discard deviant flux values. 

First the sky image is smoothed (median smoothing) 
with a 10 x 10 window. Then, for a given sky position (a, 
(5), we compare the N pixels in the data cube that have 
seen that position with the flux of the smoothed sky image. 
For most of the sky position, the N pixels are distributed 
around the smoothed sky estimate. On the other hand, 
at point sources positions, the N pixels will be generally 
above the smoothed sky estimate. Furthermore, it is also 
possible to find sky positions where most of the N pixels 
have fluxes under the smoothed sky level. Here is how we 
deal with each case: 

1. Most of the N pixels are around the smoothed sky level 
(SML). 

The new sky flux is the average of {SAIL — 3 x noise < 
I < SML + 3 x noise) 

2. Most of the N pixels are above the smoothed sky level. 
The new sky flux is the average of / > SML 

3. Most of the N pixels are under the smoothed sky level. 
The new sky flux is the average of (SML — 3 x noise < 
I < SML). 

This method allows to globally reduce the noise level and 
keep a good photometry of point sources. This is the final 
data processing step. The sky images obtained after bad 
pixels removal of the first and second GRB observations 
are presented in Figure |^ D and Figure || E. 

6. Assesment of the method 

6.1. Comparison of the two GRB observations 

From the comparison of the final maps of the two GRB 
observations (Figure || D and Figure || E) we can estimate 
the reliability of our processing. At first glance we see 
that the structure of the diffuse emission is very similar in 
both maps; the LTT correction applied seems to restore 
properly the large scale structure. Furthermore, almost 
all point like structures are present in both maps giving 
confidence in our bad pixels identification. 

The difference of the two final sky images is shown in 
Figure ^|F. It is dominated by small scale structure noise 
but large scale structures are also apparent. These are 
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memory effects are not fulfy corrected on point sources 
and that we are undersampling the PSF. One also notices 
that the noise level is higher at the edges of the differ- 
ence map, due to less redundancy in these regions. The 
standard deviation of the difference map (Figure ||F) in 
the central part of the field is 0.06 ADU/g/s. Therefore, 
the noise on each GRB final maps can be estimated at 
0.06/V2 = 0.04 ADU/g/s. 

It is clear from the difference map (Figure [sjF) and 
from the impact of the LTT on the sky image that noise 
is present at all scales. The processing presented in this 
paper affects the signal at various scales. To characterize 
the noise as a function of angular scale we use the structure 
function of second order: 



B 2 (l) 



Ulskyjr) - I sky (r + l)f 
N(l) 



(12) 



where the sum is over all the N(l) pairs of points separated 
by a distance I. To estimate the noise at each scale in the 
first map, the one affected by all instrumental effects, we 
compute the structure function on the difference between 
the first map of the GRB1 observation and the final one of 
the GRB2 observation. This difference map, where the sky 
is removed, is dominated by the noise of the GRB1 obser- 
vation. The structure function of this difference map rises 
strongly from small to large scale (see Figure Kjl), mainly 
due to the presence of the LTT. To estimate the noise at 
the end of the processing, we have computed the structure 
function on the difference between the final maps of the 
two GRB observations (see figure ^F) . This time the struc- 
ture function is very flat (see figure ^) . The noise level is 
reduced by a factor ten at 8 arcmin scale and a factor two 
at the resolution limit. This indicates that the noise level 
has been lowered at all scales and that it is now uniform 
at all scales. 

6. 2. Study of the noise sources 

The goal of this section is to show that the high spatial 
frequency noise of our maps is close to the optimal value 
obtained with stabilized ground calibration data with no 
glitches. The data are affected by many sources of noise. 
First, there are the classical quantum photon noise and 
the detector readout noise which have been extensively 
studied in the pre-launch calibration phase (Perault et al. 
1994). A conservative value of the readout noise is given 
in the ISOCAM cookbook: 1.5 ADU/g. Secondly, memory 
effects (short term transient, long term transient and slow 
glitches) and fast glitches with small amplitude increase 
substantially the noise levels. These non-gaussian events 
may prevent to reach the optimal sensitivity. 

The noise level is measured on the flux history of pix- 
els for which fast glitches have been removed. We have 
selected pixels not affected by slow glitches. We quantify 
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Fig. 9. Structure function of the difference between the 
two GRB final maps (bottom line) and of the difference 
between the first map of GRB1 and the final map of GRB 2 
(top line). 

Table || that the noise level of the two GRB observations 
is in total agreement with the readout and photon noise 
estimated from calibration data (for a 41.5 ADU/g/s flux). 

The flux I s k v at a given position in the final sky image 
is the average of N independant flux measurements. The 
error 8 on I s k y can be estimated by 

(13) 



N 



where a is the standard deviation of the N flux measure- 
ments used to compute I s k y ■ We have computed the error 
map at each step of the processing. In Figure 10 is shown 
the error 5 of each sky position for the sky image obtained 
before the LTT correction (Figure |To|A) and for the final 
sky image (Figure [lO|B) of the GRB1 observation. In the 
Figure [TojA one sees an enhancement of the error in the 
southern part of the image, due to the presence of the 
LTT. Glitches and periodic flat defects appear as noise 
peaks in this error map. The structure of the final error 
map (Figure fl"o| ) is dominated by the redundancy effect: 
the noise is higher on the edges of the map as less flux 
measurements were obtained in these regions. 

For each sky image of Figure ||, Table || lists the median 
error S, the median redundancy N and the median stan- 
dard deviation of the N flux measurements averaged at 
each sky position. The first thing to notice is that the noise 
level decreases gradually through the processing. In the fi- 
nal sky image, the median a is only 5% above the noise 
calibration measurement for the GRB1 observation. The 
noise in the GRB2 observation is exactly the one obtained 
with stabilized ground calibration data with no glitches. 

The dispersion of the difference between the tw o fi- 
nal maps (divided by y/2) is 0.04 ADU/g/s (see § O 
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Table 2. Noise table. The upper panel gives the noise level computed on the flux history of a pixel. The bottom one 
indicates the noise level in the sky image at each step of the processing. 





Noise on a single ISOCAM readout (readout and photon 


noise] 




Ground alibration noise estimate (for / = 41.5 ADU/g/s) 






1.15 


Noise measured on the flux history of a pixel (GRB1) 






1.15 




Sky Image 


Figure 


S a 


N b 


a 








(ADU/g/s) 


O-fy) 




(ADU) 


GRB1 


Before LTT correction 


8A 


0.062 


15.0 


76 


2.72 c 




After LTT correction 


8B 


0.053 


15.8 


76 


2.32 c 




Variable Flat Field 


8C 


0.046 


11.13 


76 


2.01 c 




Final map 


8D 


0.031 


7.50 


58 


1.21 c 


GRB2 


Final map 


8E 


0.030 


7.26 


60 


1.15" 



a median error value of each sky image. 
b median redundancy of each sky image. 



c a = Error x ^/Redundancy x Integration time. 



is partly due to the increase noise on point sources and 
to the imperfection of the LTT correction. Nevertheless, 
considering the amplitude of the instrumental effects (3 
ADU/g/s for the LTT) we think that the comparison be- 
tween the two GRB observations is a strong validation of 
the whole processing. 

Finally we conclude from the numbers of Table || that 
the noise level in our final maps are dominated by the 
readout and photon noise. The other instrumental effects 
are corrected in the data reduction. Such ISOCAM sensi- 
tivity has been obtained in recent point sources extraction 
studies (Desert et al. 1999; Aussel et al. 1999) where the 
low frequency diffuse emission is removed. It is the first 
time that such a sensitivity is reached for the emission 
structure on all scales. 

7. Conclusion 

The data processing techniques presented here are based 
on the fact that each position on the sky has been observed 
several times and by different pixels along the observation. 
Taking advantage of this redundancy information inherent 
to raster mode observations, we are able to control some 
instrumental effects which degrade the signal, reduce the 
noise level and finally bring images obtained with the LW 
channel of ISOCAM to a quality level close to optimal. 
The principal of the method presented could be general- 
ized to all observations with spatial redundancy. 

The quality of the final images is limited by our abil- 
ity at separating the various instrumental effects. We have 
adopted a sequencial approach where each of the problems 
is treated one after the other, but a global minimization 
would be a better way to find the optimal solution. How- 
ever such an approach seems to be out of reach for ISO- 
CAM since several instrumental problems (e.g. the LTT 
and response variation after glitches and point sources) 
are presently not understood. 



solar system extended objects (comet dust trails), nearby 
clouds and star forming regions, images from diffuse emis- 
sion in the Galactic plane and external galaxies. Several 
publications based on the interpretation of ISOCAM ob- 
servations reduced by the present method are being pre- 
pared. 
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Appendix A: Flat fielding 

A.l. Averaging on a sliding window 

One way to estimate a flat field for each readout of a data cube 
consists to compute a sliding window average of the data cube: 



F(x,y,U) = 



ti+N/2 

- y 

tA=ti-N/2 



Iobs(x,y,tj), 



(A.l) 



where N is the size of the sliding window. 

The variable flat field is normalized to 1 over the 11x11 
central part of F(x,y,i): 



F(x,y,U) 



F(x,y,U) 



{F(10 : 21,10 : 21,£ i )>' 



(A.2) 



Here we suppose that during the N readouts the camera 
has observed many positions on the sky and each pixel of the 
detector has seen the same flux in average. This is a good 
approximation if the sky observed is rather smooth and does 
not contain too much small scale structures. If the emission is 
fairly uniform on the detector, the number of images required 
to comDute the flat field mav be relatively small. The choice of 
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use N = 100). Extreme high and low fluxes pixels were re- 
jected (top and bottom 15%) in order to exclude glitches and 
point sources. The main limitation of this method is that sky 
structures may be still present in the variable flat field. 

A. 2. Perturbated single flat field 

The perturbated single flat field method is based on the follow- 
ing assumption. We suppose that the variable flat field is, to 
first order, a single flat field F(x,y) computed on the whole 
data cube (see § The temporal variations of the flat field 
are treated as perturbations of the single flat, and we have: 



I obs (x,y,t) = [1 + S(x,y,t)]F(x,y)I s k v (x,y,t), 



(A.3) 



where 8(x, y, i) is the perturbation term that take into account 
flat field deformations. 

Flat deformations are mainly due to slow glitches and bad 
short term transient correction (especially on point sources). 
Therefore, for one given readout, they are considered as small 
scale defects, and 5{x,y,i) is dominated by high frequency 
structures in the (x, y) space. 

On large scale, the quantity 



= Isky(l + S{x,y,t)) 



(A.4) 



I obs (x,y,t) 
F(x,y) 

is dominated by real large scale structures (I s k y (x,y,t)) and 
not by flat defects (S(x,y,t)). 

The low frequency sky emission IsLF(x,y,t) is esti- 
mated by smoothing (median filtering) each readout of 
Iobs(x,y,t)/Flat(x,y). The size of the smoothing window has 
to be smaller than the smallest sky structures and larger than 
the largest flat defects. In most cases a compromise must be 
found (typically a 7x7 window). Then the flat field deforma- 
tions are estimated using: 



(l + 6(x,y,t)) 



Iob s (x,y,t) 



F{x,y)IsLF(x,y,t)' 



(A.5) 



We get rid of residual high frequency rea l sky structures by 
applying a temporal sliding ave rage (see § A.l) over the right 
hand side term of equation A.5. 



The variable flat field is finally obtained with: 
Flat(x,y,t) = [l + 5(x,y,t)]F(x,y) 



(A.6) 



The averaging on a sliding window method (see § A.l) follows 



low frequency temporal flat field deformations but, as a limited 
number of readouts are averaged together, sky structures are 
still present in the flat field. The perturbated single flat field 
is a better estimate since low frequencies of the incident flux 
cube are removed. 



We have evaluated the statistical error on the offset correc- 
tion A(ij): 



°{U) = J n{ n_ I) £ (A(*0-*(M,fc,*i)) 2 (B.2) 

where n is the number of pixel pairs. 
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Appendix B: LTT error determination 

The error on the LTT correction is estimated using the stan- 
dard deviation constructed with all the pixel pairs used to com- 
pute A(t). The offset correction is based on Equation^. For a 
pixel pair ((a,S,ti) and (a, 5, tj)) that have seen the same sky 
at different time, the offset at time i, is given by: 
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Fig. 8. LW10 images of a the first GRB970402 observation after deglitching, dark substraction, transient correction 
(A), long term transient correction (B), variable flat field (C) and bad pixels removal (D). Image (E) is the final map 
of the second GRB970402 observation and image (F) is the difference between (D) and (E). For these observations, 1 
ADU/g/s corresponds to 0.242 mJy/pix or 0.286 MJy/sr. 
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Fig. 10. Noise maps before (left) and after (right) the processing. For this observation 1 ADU/g/s corresponds to 
0.242 mJy/pix. 



